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ABSTRACT 


In recent years our understanding of the dense matter equation of state (EOS) of neutron stars has 
significantly improved by analyzing multimessenger data from radio/X-ray pulsars, gravitational wave 
events, and from nuclear physics constraints. Here we study the additional impact on the EOS from 
the jointly estimated mass and radius of PSR J0740+6620, presented in Riley et al. (2021) by analyzing 
a combined dataset from X-ray telescopes NICER and XMM-Newton. We employ two different high- 
density EOS parameterizations: a piecewise-polytropic (PP) model and a model based on the speed of 
sound in a neutron star (CS). At nuclear densities these are connected to microscopic calculations of 
neutron matter based on chiral effective field theory interactions. In addition to the new NICER data 
for this heavy neutron star, we separately study constraints from the radio timing mass measurement 
of PSR J0740+6620, the gravitational wave events of binary neutron stars GW190425 and GW170817, 
and for the latter the associated kilonova AT2017gfo. By combining all these, and the NICER mass- 
radius estimate of PSR J0030--0451, we find the radius of a 1.4 Mc; neutron star to be constrained 
to the 95% credible ranges 12.33* 0.76 km (PP model) and 12.18*028 km (CS model). In addition, we 
explore different chiral effective field theory calculations and show that the new NICER results provide 
tight constraints for the pressure of neutron star matter at around twice saturation density, which 
shows the power of these observations to constrain dense matter interactions at intermediate densities. 


Keywords: dense matter — equation of state — stars: neutron — X-rays: stars — gravitational waves 


1. INTRODUCTION 


Our understanding of the dense matter equation 
of state (EOS) of neutron stars has made significant 
progress over the last few years due to the arrival of 
new avenues to measure observables like mass, radius 
and tidal deformability, that connect to the behavior 
of matter at supranuclear densities. Recently NASA's 
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X-ray timing telescope, the Neutron Star Interior Com- 
position Explorer (NICER), has delivered the first joint 
measurement of mass and radius through pulse profile 
modeling of the millisecond pulsar PSR J0030+0451 
(Riley et al. 2019; Miller et al. 2019). The impact of 
this measurement on the dense matter EOS has been 
extensively studied in various EOS frameworks (see, 
e.g., Raaijmakers et al. 2019; Miller et al. 2019; Raai- 
jmakers et al. 2020; Essick et al. 2020b; Landry et al. 
2020; Dietrich et al. 2020; Jiang et al. 2020; Al-Mamun 
et al. 2021), including EOS with phase transitions to 
quark matter (see, e.g., Xie & Li 2021; Li et al. 2020; 
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Tang et al. 2021; Blaschke et al. 2020; Alvarez-Castillo 
et al. 2020) and models that explore the possibility of 
there being two stable neutron star branches (Christian 
& Schaffner-Bielich 2020). 

Concurrently, the second and third observing runs of 
LIGO/Virgo have so far resulted in the confirmed gravi- 
tational wave detections of two (most-likely) binary neu- 
tron star mergers: GW170817 (Abbott et al. 2017c, 
2019a) and GW190425 (Abbott et al. 2020a). By ac- 
curately measuring the gravitational wave phase, limits 
can be put on the EOS-dependent tidal deformability of 
the neutron stars (Flanagan & Hinderer 2008; Hinderer 
et al. 2010). While for GW170817 the tidal deformabil- 
ity could be measured within a 90% highest posterior 
density interval when adopting low spin priors (see, e.g. 
Abbott et al. 2018, 2019b), the low signal-to-noise ratio 
(SNR) of GW190425 resulted in only weak upper limits 
on the tidal deformability even when assuming low spins 
(Abbott et al. 2020a). We consider the ~ 2.6 Mọ sec- 
ondary object in GW190814 (Abbott et al. 2020d) to be 
a black hole (Nathanail et al. 2021), and will therefore 
not use this third event in our analysis. 

At nuclear densities, the EOS is well constrained by 
nuclear theory and experiments (see, e.g., Tsang et al. 
2012; Lattimer & Lim 2013; Huth et al. 2021). In par- 
ticular, many-body calculations based on chiral effective 
field theory (EFT) interactions have enabled system- 
atic predictions for the neutron matter EOS up to nu- 
clear saturation density including theoretical uncertain- 
ties (see, e.g., Hebeler et al. 2013; Tews et al. 2013; Lynn 
et al. 2016; Drischler et al. 2019; Drischler et al. 2020). 
Up to saturation density, the resulting symmetry energy 
and pressure of neutron matter are also consistent with 
extractions from nuclear experiments (Lattimer & Lim 
2013), including from measurements of the dipole polar- 
izability of neutron-rich nuclei (Roca-Maza et al. 2015; 
Birkhan et al. 2017; Kaufmann et al. 2020). Taking 
these results at nuclear densities, combined with stan- 
dard crust EOS, different extrapolations to high densi- 
ties have been found to lead to NS radii consistent with 
all multimessenger observations (see, e.g., Raaijmakers 
et al. 2020; Essick et al. 2020b; Annala et al. 2020; Diet- 
rich et al. 2020; Biswas et al. 2021). Recently, the results 
of PREX-II have pointed to higher pressures (Adhikari 
et al. 2021; Reed et al. 2021), but with very large uncer- 
tainties, so that in a combined analysis with astrophys- 
ical and chiral EFT constraints, the overall consistency 
still persists (Essick et al. 2021). 

NICER data has now enabled a joint estimate 
of the mass and radius of the high-mass rotation- 
powered millisecond pulsar PSR J07404-6620. Since 
PSR J07404-6620 (unlike PSR J0030--0451) is in a bi- 


nary with an inclination that allows measurement of the 
Shapiro delay, its mass can be measured independently 
via radio timing. Cromartie et al. (2020) reported a 
mass of 2.14 7-10 Mo, and a joint campaign by the 
North American Nanohertz Observatory for Gravita- 
tional Waves (NANOGrav) and the Canadian Hydro- 
gen Intensity Mapping Experiment (CHIME)/Pulsar 
collaborations has now resulted in an updated mass of 
2.08 + 0.07 Mo (Fonseca et al. 2021). 

Riley et al. (2021) have used this mass measurement as 
an informative prior for pulse-profile modeling analysis 
that is joint over the phase-resolved spectroscopic data 
from NICER and phase-averaged data from the XMM- 
Newton European Photon Imaging Camera (EPIC). The 
inclusion of the smaller XMM-Newton (hereafter XMM) 
data set allows for better constraints on the propor- 
tion of the X-ray emission that is attributable to back- 
ground rather than PSR J0740+6620, ultimately acting 
to cut out solutions with high compactness. This re- 
sults in an inferred radius of 12.39 525 km, and a mass 
of 2.072* 0097 Mo that is little changed from the radio 
prior. For a full description of the methodology em- 
ployed in the mass-radius inference we refer the reader 
to Riley et al. (2021). 

In this Letter, we use the mass and radius from Ri- 
ley et al. (2021) for PSR J07404-6620 as input for in- 
ferring the dense matter EOS, combining it with other 
constraints from nuclear theory and multi-messenger ob- 
servations. It should be considered as a follow-up to 
our previous work that built on NICER's results for 
PSR J0030+0451 (Raaijmakers et al. 2019, 2020), where 
in this work we explore also a broader range of multi- 
messenger constraints. As the high-density constraints 
from astrophysical observations get more precise, with 
the new NICER results and future LIGO/Virgo mea- 
surements, it will be intriguing to see them play out 
with the present nuclear constraints. In this Letter, we 
also explore this for the new NICER results and how 
they constrain the EOS above nuclear densities starting 
from different chiral EFT calculations!. 


2. INFERENCE FRAMEWORK 


In this work we will closely follow the analysis frame- 
work developed previously in Greif et al. (2019), Raaij- 
makers et al. (2019) and Raaijmakers et al. (2020). Be- 
low, we summarize this method and highlight several 
updates to the framework. 


1 The posterior samples and scripts to make the plots in this 
Letter are available in à Zenodo repository at Raaijmakers et al. 
(20212). 
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We consider two EOS parameterizations: i) a piece- 
wise polytropic (PP) model with three segments be- 
tween varying transitions densities (Hebeler et al. 2013), 
and ii) a speed-of-sound (CS) model first introduced 
in Greif et al. (2019). To capture the uncertainty 
in the EOS around nuclear saturation density (no = 
0.16 ia), both parameterizations are matched to a 
power law fit of a range of EOS calculated from chiral 
effective field theory interactions (Hebeler & Schwenk 
2010; Hebeler et al. 2013) below 1.1no. At densities be- 
low 0.5no this power law fit is connected to the BPS 
crust EOS (Baym et al. 1971). 

To constrain these EOS parameterizations, governed 
by the EOS parameters 0, we employ Bayes’ theorem 
and write the posterior distributions of the EOS param- 
eters and central energy densities € as 


p(0, £ |d, M) c p(@|M) p(e|0, M) p(4|0, M), (1) 


where M denotes the model including all assumed 
physics and d the dataset used to constrain the EOS, 
consisting of, e.g., radio-, X-ray and gravitational wave 
data. When assuming each of these datasets to be inde- 
pendent of each other, we can separate the likelihoods 
and write 


p(0, €| d, M) œ p(0 | M) p(e |8, M) 


x I»^. Aoi, Mi i, Moi | daw,i(, dem,i)) 
x |] (Vj, R; | dxicur,j) 
j 


x LH [r | dradio,k) . (2) 
k 


Here the products run over the number of different ob- 
served stars, or mergers, in the case of the gravita- 
tional wave data. Furthermore, in Equation (2) we have 
equated the nuisance-marginalized likelihoods to the 
nuisance-marginalized posterior distributions derived in 
Riley et al. (2019); Fonseca et al. (2021); Riley et al. 
(2021); Abbott et al. (2019a, 2020a). This approxima- 
tion is justifiable when the priors used in estimating 
these nuisance-marginalized posterior distributions are 
uninformative, which for simplicity we will assume to 
be a uniform prior in this case. The posterior distri- 
butions derived by Riley et al. (2019) and Riley et al. 
(2021) already use a jointly uniform prior in mass and 
radius. The posterior distributions derived by Abbott 
et al. (2019a) and Abbott et al. (2020a) use a jointly 
uniform prior in the tidal deformabilities of the two 
components A; within the range A; C [0,5000] (for 
GW190425 the upper bound of As was set to 104.). 


The prior on the detector frame masses, which are red- 
shifted with respect to the source frame masses (Maec = 
M;(1+2)), is uniform within the range Mae C [0.5, 7.7] 
and Mae C [1, 5.31] for GW170817 and GW190425 re- 
spectively. However, the posterior distribution on com- 
ponent masses from gravitational waves is highly de- 
generate because of the accurately measured chirp mass 
M. = (MLM3y3/5 /(Mi + M3)!/5. To speed up the 
convergence of our parameter estimation, we therefore 
transform the gravitational wave posterior distributions 
to include the two tidal deformabilities, chirp mass and 
mass ratio q, while reweighing such that the prior distri- 
bution on these parameters is uniform. Further, we also 
fix the chirp mass to its median value, since the small 
uncertainty in this parameter does not affect the EOS 
parameter estimation (see Raaijmakers et al. 2020), and 
thus have: 


p(0, € | d, M) œ p(8] M) p(e |0, M) 


x [[eQui ^2; qi | Me, daw i(, dem,i)) 
x L[»Qu. R; | dicen) 
j 


x Whee | dradio,k) : (3) 
k 


Fixing the chirp mass means that the vector € only con- 
tains one central density per merger, where the tidal 
deformability of the second component is now set by 
Ag = A»(0;q). If a gravitational wave event has an 
associated electromagnetic (EM) counterpart, the likeli- 
hood for that event becomes a product of the nuisance- 
marginalized posterior distribution from the gravita- 
tional wave data and the nuisance-marginalized poste- 
rior distribution from the EM analysis, such that: 


p(Ai, Ae, q| Me, daw,dam) x p(Ai, ^2, q| Mc, daw) 
x pA, A2,q| Me, dem). (4) 


Obtaining the posterior distribution p(A1, A», q | Me, dem) 
is discussed in Section 3.2.1 for the specific case of 
AT2017gfo, the kilonova associated with GW170817 
(see, e.g., Abbott et al. 2017a,b; Arcavi et al. 2017; 
Coulter et al. 2017; Chornock et al. 2017; Cowperth- 
waite et al. 2017; Kasliwal et al. 2017; Nicholl et al. 
2017; Tanvir et al. 2017). 

We then sample from the posterior distribution 
p(@,e|d,M), compute the corresponding M, R, and 
A, and then evaluate the likelihood by applying a kernel 
density estimation to the posterior distributions from 
Riley et al. (2019, 2021); Abbott et al. (2019a, 2020a) 
using the nested sampling software MULTINEST. The 
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Figure 1. Constraints on the mass-radius relation of neutron stars, given the posterior distribution on EOS parameters 0 using 
the PP model (left) and CS model (right panel). The constraints from the updated radio timing mass of PSR J07404-6620 from 
Fonseca et al. (2021) (present work, green) are compared to the mass from Cromartie et al. (2020) used in our previous works 
(Raaijmakers et al. 2019, 2020) (orange, dashed-dotted), showing both the 68% and 95% credible regions. The black dashed 
lines indicate the 9596 credible region of the prior distribution. Note that the slightly lower mass measurement does not have a 


significant impact on the EOS posterior. 


same prior distribution p(@|M) is used as in previous 
work; we refer the reader to Section 2.3 of Raaijmakers 
et al. (2020) and references therein for a more detailed 
description. 


3. EOS CONSTRAINTS 


In this Section we investigate the impact of the 
Riley et al. (2021) mass-radius measurement for 
PSR J0740--6620 on the dense matter EOS, both sepa- 
rately and when combined with previous constraints. 


3.1. Radio mass measurement of PSR J0740-- 6620 


Firstly, we constrain the EOS using the updated mass 
measurement of 2.08 + 0.07 Me for PSR J0740--6620 
derived using radio timing (Fonseca et al. 2021), and 
compare this to the constraints from the previously pub- 
lished mass of 2.14*049 Mo (Cromartie et al. 2020). In 
Figure 1 we show the posterior distribution on EOS pa- 
rameters 0 when transformed to the mass-radius param- 
eter space. We note that, as expected, the slightly lower 
updated mass measurement shifts the posterior distribu- 
tions to lower maximum neutron star masses and lower 
radii, although the effect is almost negligible. Since the 
radio timing mass measurement is already incorporated 
in the joint mass-radius estimate from NICER we will 
not use this measurement in the remainder of this work. 


3.2. GW170817 and GW190425 


The gravitational wave events GW170817 and GW190425 


have so far been the only confirmed neutron star bi- 
nary mergers during the recent observing runs of the 
LIGO/Virgo collaboration (Abbott et al. 2020c). Al- 
though both events have a non-negligible chance of 
being neutron star-black hole mergers (see, e.g., Yang 
et al. (2018); Ascenzi et al. (2019); Coughlin & Di- 
etrich (2019); Hinderer et al. (2019) for GW170817 
and e.g., Kyutoku et al. (2020); Han et al. (2020) for 
GW190425), in the following we will assume both ob- 
jects to be neutron stars. We use the low-spin? posterior 
distributions on tidal deformability and mass ratio with 
the IMRPhenomPv2 NRTidal? waveform model (Hannam 
et al. 2014; Khan et al. 2016; Dietrich et al. 2019) 
for GW170817 and GW190425. Furthermore we use 
the median chirp mass values of M. = 1.186 Mo for 
GW170817* and Me = 1.44 Mo for GW190425°. 

'The upper panels of Figure 3 show the posterior dis- 
tributions on the EOS for both events in the mass-radius 
space. We note that the constraints on tidal deforma- 


? The low-spin assumption is chosen to be consistent with mea- 
surements of spins in Galactic neutron star binaries that merge 
within a Hubble time. 

3 See Table 1 of Abbott et al. (2019a) for a description of the 
waveform model. 

4 https:/ /dcc.ligo.org/LIGO-P1800370/public/ 

5 https: //dec.ligo.org/LIGO-P2000223/public/ 
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Figure 2. In blue we show the bolometric luminosity of 
GW170817 from the data compiled in Kasliwal et al. (2017). 
The red band contains 9596 of the light curves of the posterior 
distribution when fitted with the model described in Section 
3.2:1. 


bility from GW170817 give more support to softer EOS, 
although the 9596 credible region spans a relatively large 
range of radii. GW190425 only led to weak upper limits 
on the tidal deformability due to its low SNR and single- 
detector detection. The EOS is however constrained as a 
result of the high mass of the primary component (with 
95% credible range 1.60— 1.87 Mo), excluding EOS that 
do not support these masses. 


3.2.1. AT2017gfo 


Following the detection of GW170817 an EM counter- 
part was observed across the frequency spectrum (see, 
e.g., Abbott et al. (2017a,b) and references therein; 
Coulter et al. (2017); Chornock et al. (2017); Drout 
et al. (2017); Hallinan et al. (2017); Kasliwal et al. (2017, 
2019); Margutti et al. (2017); Pian et al. (2017); Smartt 
et al. (2017); Troja et al. (2017)). Of particular interest 
here is the thermal infrared-optical-ultraviolet transient 
powered by radioactive decay of r-process nucleosynthe- 
sis in the neutron-rich material ejected during merger; 
the so-called kilonova or macronova (e.g. Li & Paczyński 
1998; Kulkarni 2005; Metzger et al. 2010). The kilonova 
properties depend on the mass, velocity, and composi- 
tion of the ejected material, which in turn depend on 
the binary progenitor parameters such as the tidal de- 
formability of the neutron stars. Using this connection 
it is possible to constrain the EOS from the kilonova 
light curve (see, e.g., Coughlin et al. 2018; Radice & Dai 
2019; Hinderer et al. 2019; Capano et al. 2020; Dietrich 
et al. 2020). 


Here we analyze the bolometric luminosity of GW170817 


(as compiled in Kasliwal et al. 2017) via the new 


Bayesian framework outlined in Raaijmakers et al. 
(2021b). We consider a two-component kilonova model, 
where the first component, the dynamical ejecta, is asso- 
ciated with material ejected through tidal forces and the 
shock interface between the two neutron stars (see, e.g., 
Radice et al. 2018, and references therein). The second 
component is associated with neutrino-driven winds or 
material ejected through viscous forces. We connect the 
outflow properties of these components to the binary 
progenitor properties by using the formulae presented 
in Krüger & Foucart (2020) for dynamical ejecta and 
disk mass, which are fitted to numerical simulations of 
compact mergers. The velocity of the dynamical ejecta 
is calculated using the formula in Coughlin et al. (2019), 
while the velocity of the disk wind ejecta is left as a free 
parameter. The dynamical ejecta includes both mate- 
rial ejected through tidal forces and material ejected 
through shocks on the contact interface between the 
stars (see, e.g., Sekiguchi et al. 2016; Dietrich & Ujevic 
2017; Tanaka et al. 2020; Nedora et al. 2021). To dis- 
tinguish these we consider two different opacities in the 
dynamical ejecta, corresponding to tidal tail and shock 
ejecta , where the latter is less neutron-rich compared 
to the tidal tail and thus has a lower opacity (see Table 
1). For simplicity we take a single opacity for the disk 
wind ejecta. The outflow properties are then connected 
to a bolometric luminosity through the semi-analytic 
light curve model by Hotokezaka & Nakar (2020). The 
priors on all parameters are shown in Table 1. 

The fit to the bolometric luminosity of AT2017gfo us- 
ing the data compiled in Kasliwal et al. (2017) is shown 
in Figure 2, showing all datapoints to be contained 
within the 9596 credible region of the posterior distri- 
bution. In the lower panels of Figure 3 we show the 
updated prior distribution for the EOS with GW170817 
and with the inclusion of AT2017gfo. The EM data 
gives more posterior support to stiffer over softer EOS, 
due to the estimated ejected mass requiring a neutron 
star with larger tidal deformability. The estimated ra- 
dius of a 1.4 Mc neutron star for the PP and CS model is 
12.12* 110 km and 11.53*1 19 km, respectively, which is 
broadly consistent with multimessenger constraints ob- 
tained by other works (see, e.g., Coughlin et al. 2019; 
Dietrich et al. 2020; Capano et al. 2020; Breschi et al. 
2021; Nicholl et al. 2021). Important to note is that the 
EM modeling of the kilonova here is simplified and relies 
on a few assumptions that are known to affect results, 
such as spherical ejecta geometry (see, e.g., Heinzel et al. 
2021; Korobkin et al. 2021), fixed nuclear heating rate 
(see, e.g., Barnes et al. 2020), and an incomplete map- 
ping between properties of the binary system and the 
ejecta outflows. It is also dependent on the choice of 
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Figure 3. Upper panels: Constraints on the mass-radius relation of neutron stars, given the posterior distribution on EOS 
parameters 0 using the PP model (left) and CS model (right) when analyzing the gravitational wave events GW170817 (Abbott 
et al. 2017c) and GW190425 (Abbott et al. 2020a), both separately and combined. The estimated tidal deformability from 
GW170817 offers more posterior support for softer EOS, and thus lower radii. For GW190425 only weak upper limits could 
be set on the tidal deformability, but the relatively high estimated mass of the primary object disfavors softer EOS, as we are 
not considering any high-mass information from radio pulsars here. Lower panels: The change in the posterior distribution 
on the EOS when including information from the kilonova associated with GW170817, AT2017gfo (Kasliwal et al. 2017). The 
estimated mass that was ejected during the merger favors higher tidal deformabilities, and thus constrains the mass-radius space 


at low radii. 

light curve modeling, where the distinction can be made 
between semi-analytic modeling (such as in this work 
and, e.g., Breschi et al. 2021; Nicholl et al. 2021) and in- 
terpolating between radiative transfer simulations (e.g., 
Coughlin et al. 2018; Dietrich et al. 2020). We use a 
semi-analytical model from Hotokezaka & Nakar (2020), 
which for the current statistical uncertainty in gravi- 


tational wave parameter estimation and uncertainty in 
light curve observations produces consistent results to 
full radiative transport models (see, e.g., Coughlin et al. 
2020a,b), although this will change in the future with 
improved gravitational wave detectors and optical tele- 


scopes. 
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Figure 4. Upper panels: The 68% and 95% credible regions of the EOS given the mass-radius estimate of PSR J0740+6620 
by Riley et al. (2021), using the PP model (left) and CS model (right). The black dashed lines and orange dashed-dotted lines 
indicate the 95% credible region of the prior and the constraints given the radio mass measurement of PSR J0740+6620 by 
Fonseca et al. (2021), respectively. The red contour shows the posterior distribution on central energy density and pressure for 
this source, and in the inset we plot the KL-divergence as a function of energy density. Lower panels: Same as upper panels but 
for the mass-radius space. Also shown in blue, dotted lines is the 95% credible region of the EOS posterior distribution, when 
analyzing the result from Riley et al. (2021) without the inclusion of the XMM-dataset (so NICER only). In addition, we show 
the mass-radius posterior for PSR J0740+6620 by Riley et al. (2021) as dark-green contours (68% and 95%). Note that when 
considering both NICER and XMM data, the posterior distribution (green shaded) is very close to the constraints obtained 
from the radio mass measurement of PSR J0740+6620 (orange), due to this mass-radius posterior (dark green) showing support 
over an extended range of radii. 


3.3. NICER mass-radius and multimessenger (2021). They find a radius of 12.30" 5 35 km, and a 
constraints mass of 2.072*0067 Mo, where the upper and lower 


limit bound the 68% credible regions. The EOS re- 
sults are shown in Figure 4, both in energy density- 
pressure and mass-radius space. From the Kullback- 


Next we study the constraints on the EOS from the 
new mass-radius estimate of PSR J0740+6620 using 
data from NICER and XMM, presented in Riley et al. 
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Figure 5. Upper panels: Constraints on the mass-radius space of neutron stars, given the posterior distribution of EOS 


parameters 0 using the PP model (left) and CS model (right). Shown are the 68% and 95% credible regions when analyzing 
PSR J0030+0451, PSR J0740+6620 and the combination of the two pulsars. Note that the distribution of PSR. J0030+0451 
is different than in Raaijmakers et al. (2019), because here we have not included any high-mass pulsar information. Lower 
panels: Similar to upper panels, but when analyzing jointly mass-radius estimates from PSR J0740--6620 (Riley et al. 2021), 
PSR J0030+0451 (Riley et al. 2019), mass-tidal deformability estimates from GW170817 (Abbott et al. 2019a) and GW190425 
(Abbott et al. 2020a) and the kilonova data of Kasliwal et al. (2017) as described in Section 3.2.1. Combined, we find the radius 
of a 1.4Mo neutron star to be constrained to the 95% credible ranges 12.33*0 26 km (PP model) and 12.18*026 km (CS model). 
To show the impact of the radius measurement of PSR J0740--6620 we also plot the posterior distribution when analyzing 
combined constraints with only the 2.08Mọ mass measurement of PSR J07404-6620 (orange dashed-dotted lines). 


Leibler (KL)-divergence (Kullback & Leibler 1951) plot- 
ted as a function of energy density in the upper insets, 
we find that especially at higher energy densities there 
is a significant information gain from prior-to-posterior. 
Note that similar but, especially for the CS model, 
broader constraints are found for the posterior distri- 


bution when only using the radio mass measurement 
of PSR J0740+6620, as indicated by the orange dashed- 
dotted lines. This is a result of the mass-radius estimate 
of PSR J07404-6620 being very consistent with our prior 
ranges informed by low-density chiral EFT calculations. 
'The chiral EFT calculations do exclude however stiffer 
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Table 1. The parameters used in the model described 
in Section 3.2.1 and their prior support in the analysis of 
AT2017gfo. The notation U(a,b) here means uniformly 
drawn between boundaries a and b. 


Parameters Prior density and support 
Binary properties 

Me [Mo] U(1.18, 1.2) 

q ~ U(0.2, 1) 

Ay U(0, 2500) 

nv ~ U(0, 2500) 


Ejecta and light curve properties 


Mays [Mo] Eq. (2) Raaijmakers et al. (2021b) 
Vayn [c] Eq. (D5) Coughlin et al. (2019) 
Umin,dyn [C] ~ U(0.1, 1.0) Vayn 

Umax,dyn [c] ~ U(1.5, 2.5) Vayn 

vs [c] ~ U(vmin,dyn, Vmax,dyn) 

row (em? g^] ~ U(O.1, 5) 

Knigh [cm? g^ !] ~ U(5, 30) 

Mwina [Mo] Eq. (10) Raaijmakers et al. (2021b) 
Vwind [C] ~ U(0.03, 0.15) 

Umin,wind [C] ~ U(0.1, 1.0) vwina 

Umax,wind [c] ~ U(1.5, 2.0) vwina 

kwina [cm? g !] ~ U(0.1, 5) 


EOS with radii > 14km, where the mass-radius pos- 
terior of PSR J0740--6620 has non-negligible posterior 
support. For the CS model this effect is stronger as ad- 
ditional constraints on the speed of sound at 1.5no in 
the CS model lead to overall smaller radii than in the 
PP model (see Section 2.3 of Raaijmakers et al. 2020). 
In the mass-radius space we also plot the EOS con- 
straints given the joint NICER mass-radius estimate ex- 
cluding the XMM data. For this analysis Riley et al. 
(2021) report a value of 11.29* 52" km for the radius, 
and 2.078*0069 Mo for the mass. As this joint mass- 
radius estimate has slightly more posterior support for 
lower radii, the corresponding EOS constraints suggest 
a softening of the EOS at high densities. These results 
should be interpreted with caution however, because the 
NICER-only analysis leads to an under-prediction of the 
background (the contribution from instrumental or as- 
trophysical background to the unpulsed component of 
the pulse profile). This results in more of the unpulsed 
component being attributed to the hot regions via high 
compactness solutions. The XMM data show that a 
larger component of the unpulsed emission must come 
from true background, eliminating these high compact- 
ness solutions and increasing the inferred radius in the 


joint NICER-XMM analysis (see also Section 4.2 in Ri- 
ley et al. (2021)). 

Finally, in Figure 5 we show the constraints on the 
EOS from PSR J07404-6620, PSR J00304-0451 (first de- 
rived in Raaijmakers et al. 2019, but here no information 
on high-mass pulsars is included) and the combination of 
the two pulsars. Note that for the combined constraints, 
most of the information comes from PSR J0740--6620, 
since the 6896 credible region of the mass-radius pos- 
terior of PSR J00304-0451 covers a broad range in 
radii that are consistent with the EOS constraints from 
PSR J07404-6620. 

In the lower panels of Figure 5 we show the com- 
bined constraints on the EOS including mass-radius 
estimates from PSR J07404-6620 (Riley et al. 2021), 
PSR, J0030+0451 (Riley et al. 2019) and mass-tidal de- 
formability estimates from GW170817 (Abbott et al. 
2019a) and GW190425 (Abbott et al. 2020a) and the 
kilonova AT2017gfo (Kasliwal et al. 2017). We find that 
especially the pulsar mass-radius estimates by NICER 
favor stiffer EOS, as well as GW170817 when the as- 
sociated kilonova AT2017gfo (Kasliwal et al. 2017) is 
included. The weak constraints from GW190425 on the 
tidal deformability are also broadly consistent with the 
constraints coming from the other sources. As a com- 
parison we show the posterior distribution when com- 
bining all analyses excluding the mass-radius estimate of 
PSR J0740--6620, but with the radio mass measurement 
of Fonseca et al. (2021). We note that the additional 
radius information on PSR J0740--6620 constrains the 
softer EOS, especially for the CS model. 


4. SENSITIVITY OF POSTERIORS TO NUCLEAR 
CONSTRAINTS AT LOW DENSITIES 


To investigate the impact of the EOS constraints 
from nuclear physics we compare our analysis of 
PSR J07404-6620 using four different chiral EFT un- 
certainty bands. All bands are based on microscopic 
calculations for pure neutron matter, which are then 
extended to neutron star matter in beta-equilibrium us- 
ing the formalism discussed in Hebeler et al. (2013). In 
order to improve the description of all employed EOSs, 
we generalized the density dependence of the energy- 
density functional [see Eq. (2) in Hebeler et al. (2013)] 
by enlarging the range of the exponent y to y € [1.2, 2.5]. 

The results from Hebeler et al. (2013) formed the ba- 
sis of our previous studies (Raaijmakers et al. 2019, 
2020). The calculations for pure neutron matter were 
initially performed in Hebeler & Schwenk (2010) us- 
ing many-body perturbation theory, while the uncer- 
tainty band results mainly from variations of the cou- 
plings involved in three-nucleon interactions. Second, in 
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Figure 6. Different chiral EFT bands for the pressure of 
neutron star matter at nuclear densities, n/no in units of 
saturation density no — 0.16fm ?, and their matching to 
the BPS crust EOS at 0.5no. The different bands are based 
on microscopic calculations of neutron matter from Hebeler 
et al. (2013), Tews et al. (2013), Lynn et al. (2016) and 
Drischler et al. (2019) and include beta equilibrium (with 
protons and electrons) following the construction in Hebeler 
et al. (2013). The four chiral EFT calculations are considered 
between 0.5no and 1.1no in the analyses presented in Section 
4. Also shown are examples of the fit we use to approximate 
the EOS within these uncertainty bands, see Eq. (5), and 
connect to the BPS crust EOS. For a comparison of the chiral 
EFT bands in pure neutron matter, see Figure 1 in Huth 
et al. (2021). 


Tews et al. (2013) the calculations for neutron matter 
were improved by including for the first time all two-, 
three-, and four-neutron interactions to next-to-next- 
to-next-to-leading order (N?LO), which are predicted 
in a parameter-free way for neutron matter (see, e.g., 
Hebeler et al. 2015; Hebeler 2021 for reviews). Third, 
in Drischler et al. (2019) the calculations were further 
optimized by improving the treatment of three-nucleon 
interactions and extending the many-body expansion to 
higher orders. In addition, the EOS uncertainty bands 
also include effects from variations of regulator scales 
in state-of-the-art nucleon-nucleon and three-nucleon in- 
teractions. In this work, we use the combined 450 MeV 
and 500MeV N?LO uncertainty bands from Drischler 
et al. (2019). Finally, we include results of Lynn et al. 
(2016). These were obtained by nonperturbative quan- 
tum Monte-Carlo simulations of neutron matter at next- 
to-next-to-leading order (N?LO). This represents a com- 
pletely different many-body method than those used for 
the other three bands, and the results of Lynn et al. 
(2016) are also based on a different set of local two- and 
three-nucleon interactions derived from chiral EFT. 


Similar to Raaijmakers et al. (2020) we approximate 
the EOS within these bands with a single polytrope P — 
Nn". However, to obtain a better fit to the additional 
bands considered here, we vary the polytropic index I 
as a function of the normalization N, 


(N I Nin) 
(Nmax — Nin) 


T(N) = (Trax Thin) H Ins (5) 
where Nmin/max and I'i,/,4* are determined by fitting 
a polytrope to the lower and upper bound of the band. 
In Figure 6 we show the four different bands for the 
pressure of neutron star matter with an example of the 
fit through each band. This shows the consistency of 
these different chiral EFT calculations, with different 
methods, interactions, and approximations. The first 
point of the band where n/no > 0.5 is matched to the 
BPS crust EOS at 0.59 via a linear interpolation. 

We study the dependence of the EOS constraints on 
the different chiral EFT bands by inferring the EOS 
from the mass-radius estimate of PSR J07404-6620 us- 
ing each band and both high-density parameterizations. 
The results are shown in Figure 7. We also show the 
9596 credible region of the updated prior distribution 
when directly joining the PP or CS high-density param- 
eterization to the crust EOS at 0.5p,,. As expected 
the chiral EFT calculations mostly exclude stiffer EOS. 
While the different chiral EFT bands yield very good 
agreement on the upper bound of the radius estimates, 
the lower bound on the radius does slightly depend on 
the chiral EFT band used, especially at lower neutron 
star masses, depending on how soft the chiral EFT band 
is (see Figure 6). 

In the lower panels of Figure 7 we also show the poste- 
rior distributions on the pressure at densities n = 1.5no 
and n = 2no above the chiral EFT bands. These 
results demonstrate that the PSR J0740--6620 mass- 
radius measurement systematically prefers higher pres- 
sures at these densities compared to the correspond- 
ing prior distributions of each chiral EFT band. Fur- 
thermore, the posteriors at n = 2n agree very well 
for all chiral EFT bands and are peaked around P ~ 
109^5dyn/cm? ~ 20 MeV/fm?. 


5. DISCUSSION 


In this Letter, we have investigated the constraints 
on the EOS posed by the new joint mass-radius 
estimate from NICER x XMM data (Riley et al. 
2021), and compared and combined with multimes- 
senger EOS constraints from radio timing, gravitational 
wave mergers and their counterparts, and the previous 
PSR J00304-0451 mass-radius estimate by NICER. In 
Table 2 we summarize the results obtained in Sections 
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Figure 7. Upper panels: 9596 credible region for the mass-radius space given the mass-radius estimate of PSR. J07404-6620 by 
Riley et al. (2021), using the PP model (left) and CS model (right). The different results correspond to using the four different 
chiral EFT calculations between 0.5 and 1.1mo as shown in Fig. 6. Moroever, the red, dashed lines correspond to the 95% 
credible region, if the PP or CS parameterization is used down to 0.5no, i.e., immediately following the BPS crust, so that no 
information from chiral EFT is used. Lower panels: Marginalized posterior distributions for the pressure P above saturation 
density, at density n — 1.5no (left) and n — 2no (right) above the chiral EFT bands. 


3 and 4 for the constraints on the radius of a 1.4,1.6 
and 1.8 Mg neutron star, as well as AR = Rə — R4, 
and the maximum mass of a non-rotating neutron star 
Mov, as well as the constraints on the central energy 
density and pressure for PSR J0740--6620. 


5.1. Implications for nuclear physics 


We have studied the sensitivity of the EOS constraints 
from PSR J0740--6620 using four different low-density 
EOS calculations from chiral EFT (see Section 4). From 


the results presented in Figure 7 and Table 2 we con- 
clude that the constraints on the EOS are only weakly 
dependent on the choice of low-density calculations, al- 
though small differences exist at lower radii. Assuming 
all four low-density calculations to be equally probable, 
we can compute the Bayes’ factor K by taking the ra- 
tio of the evidence of each MULTINEST run, and assess 
whether one model is preferred over another by the data 
of PSR J0740+6620. We list the Bayes’ factors in Table 
2, where each model is compared to using the chiral EFT 
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Table 2. Key quantities from the posterior distributions obtained in Sections 3 and 4: The radius of a 1.4Mo, 1.6Mo, and 
1.8 Mo neutron star, as well as AR = Rə — Ri.4, and the maximum mass of a non-rotating neutron star Mrov. For the analyses 
of Section 4, we also show the inferred central energy density £c, the corresponding central pressure Pe, and the Bayes’ factor 
K comparing with the model using the chiral EFT band from Hebeler et al. (2013). The first four column results are for the 
different chiral EFT bands from Hebeler et al. (2013) (Heb 13), Tews et al. (2013) (Tews 13), Lynn et al. (2016) (Lynn 16), 
and Drischler et al. (2019) (Dri 19), while all other results are for the baseline inference using Heb 13. The column *Combined 
with" refers to the NICER x XMM analysis of PSR J0740+6620, the NICER analysis of PSR J0030--0451 and multimessenger 
constraints combined, while in the column “Combined without" the NICER x XMM analysis of PSR J07404-6620 is replaced 
with just the radio mass measurement by Fonseca et al. (2021). The radii are given in km, Mrov in Mo, and ec and P. in 
g/ cm? and dyn/ cm?, respectively. The upper and lower values correspond to the 9596 credible interval. 
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K 1.00 


band from Hebeler et al. (2013). All values are close to 
one, indicating that there is no substantial support for 
one model over the other, based on the mass-radius esti- 
mate of PSR J0740+6620. These results are consistent 
with the observation that predictions for pure neutron 
matter are well constrained by modern nuclear forces 
derived within chiral EFT (Huth et al. 2021; Hebeler 
2021). 

Also shown in Table 2 are the values of AR = Rə — 
Rı.4, the difference in radius of a 2 Mọ and 1.4 Mo neu- 
tron star. As pointed out by Drischler et al. (2021), the 
value of AR, if positive, can give an indication that pos- 
sibly unusual stiffening happens at high densities. We 
find however all values to be consistent with the mean 
AR being negative, but due to the broad uncertainty no 
conclusive statements can be made. 


5.2. Implications for maximum mass 


An important quantity relating to the EOS is the max- 
imum stable mass of a non-rotating neutron star, Mroy. 


Accurate knowledge of Mroy can aid in classifying com- 
pact mergers and merger remnants. In Figure 8 we show 
posterior distributions on Mroy when analyzing the up- 
dated radio mass measurement of PSR J0740--6620, 
the joint mass-radius estimate of PSR J07404-6620 
and combining GW170817, GW190425, AT2017gfo, 
PSR J0740+6620 and PSR J00304-0451. The latter re- 
sults in 95% credible ranges for Mrov = 2.234358 Mo 
and Mrov = els a ter Mo for the PP and CS model, 
respectively. This is in agreement with values previous 
found (see, e.g., Nathanail et al. 2021, and references 
therein) when assuming the secondary component in 
GW190814 was a black hole (Abbott et al. 2020d). Note 
that the higher end of the distribution in Figure 8 is very 
dependent on our choice of parameterization, as no in- 
formation is included from sources with masses above 
2.08 Mo. One could use information on the merger rem- 
nant of GW170817 to put an upper bound on Mrov 
(see, e.g., Margalit & Metzger 2017; Shibata et al. 2017; 
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Ruiz et al. 2018), but that is beyond the scope of this 
Letter. The lower end of the distribution on the other 
hand is strongly correlated with the radio mass mea- 
surement of PSR J0740+6620. The recently lowered 
mass distribution presented in Fonseca et al. (2021) re- 
sults in slightly lower values for Mrovy compared to the 
distributions found in Raaijmakers et al. (2020). 


5.3. Systematic uncertainties and framework 
comparisons 


'The analysis presented in this Letter is conditional on 
both the modeling choices of the dense matter EOS and 
on modeling choices within each analysis of the multi- 
messenger sources considered here. The sensitivity to 
the EOS modeling is explored here by employing two 
different high-density parameterizations and four differ- 
ent low-density chiral EFT calculations (see Section 4). 
From Table 2 we conclude that the CS model systemat- 
ically predicts lower radii, as a result of the additional 
constraints on the speed of sound that are not consid- 
ered in the PP model. The discrepancy between the two 
models increases with increasing neutron star mass, as 
high-mass stars depend more sensitively on the choice 
of high-density parameterization. The two models con- 
sidered here are however not exhaustive as many more 
high-density parameterizations exist (see, e.g., Lindblom 
2018; Capano et al. 2020; O'Boyle et al. 2020). 

Furthermore, we do not consider the impact of any 
systematic effects present in estimating the posterior 
distributions on M, R and A. For example, the un- 
certainty in modeling the hot regions in pulse profile 
modeling and the effect on the EOS has been studied in 
Raaijmakers et al. (2020) by using two different models 
to fit PSR J0030--0451, which led to slightly different 
constraints. For PSR J0740--6620, different assump- 
tions and priors lead to a higher estimated radius in 
the independent analysis of Miller et al. (2021) (see the 
extensive discussion of this issue in Section 4.4 of Riley 
et al. 2021), and we refer the reader to that paper for 
an EOS analysis using those results. 9. 

Measurements in A from gravitational wave data are 
also sensitive to choice of priors and gravitational wave- 
form models (see, e.g., Kastaun & Ohme 2019; Gamba 


6 Note however that one of the main reasons for the higher in- 
ferred radius reported by Miller et al. (2021) is that they do not 
truncate the prior on radius during the pulse profile modelling 
step, which Riley et al. (2021) do (truncating above 16 km, reflect- 
ing the lack of EOS models predicting higher radii, and thereby 
lowering the computational cost by reducing the parameter space). 
In the analysis by Miller et al. (2021) the lack of prior support for 
high radii is effectively incorporated at a later stage, in the EOS 
analysis. 


et al. 2020). Lastly, many different kilonova models exist 
(see, e.g., Dietrich et al. 2020; Nicholl et al. 2021; Breschi 
et al. 2021, for recent analyses) that derive slightly dif- 
ferent constraints on the EOS due to differences in mod- 
eling assumptions on, e.g., geometry, composition and 
the connection between binary properties and outflow 
properties. 

'The inference framework employed in this Letter was 
first discussed in Riley et al. (2018) and subsequently de- 
veloped in Greif et al. (2019); Raaijmakers et al. (2019, 
2020), which also introduced the chiral EFT constraints. 
Although an exhaustive comparison with other frame- 
works is out of the scope of this work, we will briefly 
mention similarities and differences with some com- 
monly used frameworks in the field. Firstly, we make use 
of two particular high-density EOS parameterizations. 
Besides many different existing choices in these param- 
eterizations, a completely different approach is to use 
non-parametric inference involving Gaussian Processes 
(see, e.g., Landry & Essick 2019; Essick et al. 2020a; Han 
et al. 2020), or discretely sampling a set of pre-computed 
EOS (see, e.g., Capano et al. 2020; Dietrich et al. 2020). 
Secondly, we compute likelihoods by performing kernel 
density estimation on posterior samples of neutron star 
properties such as mass, radius and tidal deformability 
(see also, e.g, Miller et al. 2019; Al-Mamun et al. 2021). 
It is also possible to directly infer EOS properties from 
the observational data, for example X-ray or gravita- 
tional wave data. For the first, Riley et al. (2018) argue 
that this approach would be computationally too ex- 
pensive, while for the latter this has been done by, e.g., 
Capano et al. (2020); Dietrich et al. (2020). A slightly 
different approach is used by Hernandez Vivanco et al. 
(2020), where the likelihood is computed by interpolat- 
ing marginalized likelihoods using machine learning. 


5.4. Summary and future prospects 


In summary, the new joint mass-radius estimate of 
PSR J0740--6620 significantly constrains the EOS. For 
the PP model the information gain is mostly a result of 
the high mass of the pulsar, as the 68% credible range of 
the radius estimate exactly encompasses our prior dis- 
tribution, informed by chiral EFT calculations, in that 
mass range. For the CS model the relatively high radius 
estimate does constrain the model at lower radii on top 
of constraints coming from the mass estimate. Com- 
bined with other current observational data from grav- 
itational waves and kilonova light curves, as well as the 
NICER mass-radius estimate of PSR J00304-0451, we 
find the 95% credible ranges 12.38* 0:2» km (PP model) 
and 12.23+9:38 km (CS model) for the radius of a 1.4 M 


neutron star. 
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Figure 8. Posterior distribution of the maximum mass of a non-rotating neutron star Mrov for the PP model (left) and 
CS model (right) when considering only the radio mass measurement of PSR J0740--6620, the joint mass-radius estimate of 
PSR J0740--6620 (NICER x XMM), and when combining NICER's results on PSR J07404-6620 and PSR J0030+0451 with 
GW170817 and GW190425, and AT2017gfo. For the latter ( Combined") we find a 95% credible range for Mrov = 2.231015 Mo 
and Mrov = 2.11*072 Mo for the PP and CS model, respectively. Also shown in pink is the radio mass measurement of 
PSR J07404-6620 from Fonseca et al. (2021), as the heaviest pulsar measured to date. 


In the near future, the detailed analysis of gravita- 
tional wave events observed during the second part of 
the third observing run of LIGO/Virgo are expected to 
be published, among them a few candidate events which, 
in an initial rapid classification, were identified as con- 
taining at least one neutron star. Any measured tidal 
deformability from these gravitational waves events will 
help constrain the EOS further. There were unfortu- 
nately no EM counterparts for the potential binary neu- 
tron star or black hole-neutron star events during this 
observing run. The fourth observing run is planned to 
start next year, with the LIGO and Virgo detectors close 
to their design sensitivity and KAGRA fully joining the 
network (Abbott et al. 2020b). At design sensitivity, 
GW170817-like signals will have signal-to-noise ratios of 
100 and enable measurements of tidal deformability with 
more than three times better accuracy (Capano et al. 
2020). Subsequent further detector improvements are 
already planned for the mid to late 2020s (Abbott et al. 
2020b), and an ongoing worldwide effort is paving the 
way for next decade’s third generation detectors. These 
will improve current measurements of tidal deformabil- 
ity by a factor of ~ 10 and observe the population of 
tens to hundreds of thousands of neutron star binaries, 
with EM counterparts detectable for a fraction of them 
(Maggiore et al. 2020; Sathyaprakash et al. 2019a,b). 

In the coming months, NICER is expected to deliver 
mass-radius measurements for three additional pulsars: 


two for which independent mass constraints exist (the 
~ 1.4Mo pulsar PSR J0437-4715 and the ~ 1.9 Mọ pul- 
sar PSR J1614-2230); and the pulsar PSR J1231-1411, 
which has no independently known constraint on the 
mass. There will be an update to the inferred mass and 
radius of PSR J00304-0451, using a larger data set, tak- 
ing into account improvements to our understanding of 
the NICER instrument response, and including XMM 
data in a joint analysis (as done for PSR J07404-6620). 
'There are also good prospects for narrowing the mass- 
radius measurements for PSR J0740--6620, using mod- 
els of the NICER background. All of these promise fur- 
ther improvements to our understanding of the dense 
matter EOS. 
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